##########################################################################################
## https://sunduanchen.github.io/Scissor/vignettes/Scissor_Tutorial.html
library(data.table)
library(optparse)
library(clusterProfiler)
library(fgsea)

##########################################################################################

option_list <- list(
    make_option(c("--input_file"), type = "character"),
    make_option(c("--pathway_path"), type = "character"),
    make_option(c("--out_path"), type = "character")
)

if(1!=1){

    input_file <- "~/20220915_gastric_multiple/dna_combinePublic/images/singleCell_MUC6/Diff/DiffGene.Pit.tsv"
    out_path <- "~/20220915_gastric_multiple/dna_combinePublic/images/singleCell_MUC6/Diff/GSEA"
    pathway_path <- "~/ref/Pathway/"

}

##########################################################################################

parseobj <- OptionParser(option_list=option_list, usage = "usage: Rscript %prog [options]")
opt <- parse_args(parseobj)
print(opt)

out_path <- opt$out_path
input_file <- opt$input_file
pathway_path <- opt$pathway_path

dir.create(out_path , recursive = T)

##########################################################################################
# 读入表达矩阵
exprs <- read.table(input_file, header = T, row.names = 1, sep = "\t")

##########################################################################################

for( gmt_file in grep( "gmt" , list.files(pathway_path) , value = T ) ){

    # 读入基因注释信息
    gobp_gmt <- read.gmt(paste0( pathway_path , "/" , gmt_file))
    # 对表达数据进行差异表达分析
    ## 这里去掉了基因集前缀
    gobp_gmt$term <- gsub('GOBP_','',gobp_gmt$term)
    gobp_gmt$term <- gsub('KEGG_','',gobp_gmt$term)
    gobp_gmt$term <- gsub('WP_','',gobp_gmt$term)
    gobp_gmt$term <- gsub('REACTOME_','',gobp_gmt$term)
    gobp_gmt$term <- gsub('HALLMARK_','',gobp_gmt$term)
    gobp_gmt.list <- gobp_gmt %>% split(.$term) %>% lapply( "[[", 2)

    ##########################################################################################

    geneList <- exprs[order(exprs$log_FC , decreasing=T),"log_FC"]
    names(geneList) <- rownames(exprs[order(exprs$log_FC , decreasing=T),])

    Agsea_res <- fgsea( pathways = gobp_gmt.list, 
                       stats = geneList,
                       minSize=5,
                       maxSize=500,
                       nperm=1000)
    g2 <- Agsea_res[Agsea_res$pval < 0.05,]
    g2 <- g2[order(g2$NES,decreasing = T),]

    out_file <- paste0(out_path,"/",gmt_file,".Diff.fgsea.tsv")
    write.table( data.frame(g2)[,1:7], out_file , row.names = F , sep = "\t" , quote = F )

    Up <- g2[ES > 0][head(order(pval), n=10), pathway]
    Up <- Up[!(Up %in% c(
        "ARRHYTHMOGENIC_RIGHT_VENTRICULAR_CARDIOMYOPATHY_ARVC" , "DILATED_CARDIOMYOPATHY" , "HYPERTROPHIC_CARDIOMYOPATHY_HCM" ,
        "REGULATION_OF_NEURON_PROJECTION_DEVELOPMENT"
        ))]
    Down <- g2[ES < 0][head(order(pval), n=10), pathway]
    Down <- Down[!(Down %in% c("CARDIAC_MUSCLE_CONTRACTION" , "ALZHEIMERS_DISEASE" , "PARKINSONS_DISEASE" , "HUNTINGTONS_DISEASE" , "SYSTEMIC_LUPUS_ERYTHEMATOSUS"))]
    topPathways <- c(Up, Down)
    #tmp <- topPathways[2]
    #topPathways[2] <- topPathways[5]
    #topPathways[5] <- tmp

    out_file <- paste0(out_path,"/",gmt_file,".Diff.fgsea.pdf")
    pdf(out_file , width = 18 , height = 7)
    plotGseaTable(gobp_gmt.list[topPathways], geneList, Agsea_res, 
                  gseaParam = 0.5)

    dev.off()

}